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Abstract 

While thermostated time evolutions stand on firm grounds and are widely used 
in classical molecular dynamics (MD) simulations [1], similar methods for quan- 
tum MD schemes are still lacking. In the special case of a quantum particle in a 
harmonic potential, it has been shown that the framework of coherent states per- 
mits to set up equations of motion for an isothermal quantum dynamics [2]. In the 
present article, these results are generalized to indistinguishable quantum particles. 
We investigate the consequences of the (anti-)symmetry of the many-particle wave- 
function which leads to quantum entangled distribution functions. The resulting 
isothermal equations of motion for bosons and fermions contain new terms which 
cause Bose-attraction and Pauli-blocking. Questions of ergodicity are discussed for 
different coupling schemes. 
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1 Introduction and summary 

Classical MD simulations are performed by solving Hamilton's equations of 
motion numerically. Therefore, the internal energy of the system is conserved 
during time evolution. If the ergodic hypothesis is satisfied, a time average 
of a macroscopic observable is, under equilibrium conditions, the same as 
a microcanonical ensemble average. Hence, isoenergetic molecular dynamics 
entails a method for the calculation of microcanonical ensemble averages. 
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In the canonical ensemble, however, the internal energy of the system is not 
constant, but free to fluctuate by thermal contact to an external heat bath. 
In order to adapt classical MD simulations to the problem of the calculation 
of canonical ensemble averages, i. e. to pass over from an isoenergetic to an 
isothermal time evolution, powerful methods have been developed since the 
1980s, and they are commonly used nowadays [3,4]. Some of them, like the 
Nose-Hoover chain technique [5] or the Kusnezov-Bulgac-Bauer thermostat 
[6], are based on deterministic time- reversal equations of motion leading to 
trajectories that fill the phase space of the system according to the canonical 
thermal weight. This allows the calculation of canonical ensemble averages by 
means of molecular dynamics. 

In quantum mechanics, the problem is more involved, and the pioneering ap- 
proaches of Grilli and Tosatti [7] and Kusnezov [8] were not applied very much. 
Moreover, the quantum mechanical time evolution itself is a hard computa- 
tional problem. However, a number of approximate quantum MD methods are 
available [9] . It is an open question whether they can be modified in a manner 
appropriate to permit the calculation of quantum canonical averages. 

In a different methodological approach to isothermal quantum dynamics, we 
have shown that in the special case of a quantum particle in a harmonic oscil- 
lator potential, the framework of coherent states permits to set up equations 
of motion for an isothermal quantum dynamics [2] . Following the approaches 
of Nose and Hoover or the similar KBB-technique, time-dependent pseudofric- 
tion terms are added to the equations of motion for the parameters of coherent 
states. The dynamics of the pseudofriction coefficients is designed such that 
the desired thermal weight function is a stationary solution of a generalized 
Liouville equation on a mixed quantum-classical phase space. 

The principle of indistiguishability of identical particles adds an important 
quantum aspect to the problem of isothermal quantum dynamics. The present 
article focusses on the modifications of the method presented in [2] for a sys- 
tem of two non-interacting identical particles. The (anti-)symmetry of the 
wavefunction leads to an entangled distribution function on quantum phase 
space. Surprisingly, we find that the originally classical Nose-Hoover approach 
succeeds to allow for this quantum feature. The entanglement leads to addi- 
tional terms in the equations of motion of the pseudofricional coefficients that 
act effectively like an attractive (for bosons) or repulsive force (for fermions). 
We examine whether the modified equations of motion lead to an ergodic 
time evolution, and it turns out that the additional terms improve the overall 
ergodicity of the various schemes. In addition, standard techniques to gen- 
erate ergodicity known from the classical methods are employed and tested. 
As a general result, ergodic time evolution is achieved even for low temper- 
atures without serious difficulties. However, it is indispensable to thermalize 
both particles since the effective forces do not suffice to yield an ergodic time 
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evolution of the whole system if only one particle is coupled to a thermostat. 



2 Method and setup 



2. 1 Two-particle distribution function in a harmonic oscillator potential 



The idea of the single-particle quantum Nose-Hoover method [2] is to modify 
the equations of motion of the coherent states [10] parameters r and p (which 
are combined to the complex parameter a = t -\- ^^mhw ^^^^ thut the 
quantum weight function of a single particle in a harmonic oscillator potential, 
denoted by Wqm{r,p), is sampled in time. In precise analogy to this case, we 
determine a thermal weight function w^^\ai,a2) that permits to determine 
canonical ensemble averages of two identical particles [11]. The values of £, -|- 
and — , refer to bosons and fermions, respectively. 

We write 




a: 



for the (anti-)symmetrized two-particle wavef unctions. denotes the anti- 
symmetrizing projector (i. e., | A_ ) is a two-particle Slater determinant), S+ 
the symmetrizing projector. 

The starting point for the calculation of w^^^ is the following expression for 
the calculation of a thermal expectation value {{B)) as a phase space integral. 



((5)) = itr(Se-^5) 



Ze 
1 



( «!, q;2 I o^^e ^ \ai,a2) , 



with Z^'^^ — tr ^~ j being the respective partition function. Note that we 
have dropped the projector acting upon the ket using the idempotency of 

projectors. Decomposition of e = e ^-^"^ e ^^^"^ and a cyclic shift of the 
operators under the trace enables further simplification, taking advantage of 
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the specific properties of coherent states [12]. The following expression 



Ze JJ " » ^ ■ 



'Ae\B\A,) 

(3) 



(r,-. {Ae\B\Ae) 

finally defines w\ ' as the thermal weight of the expectation value i^aJIas) • 

wf^{ai, is not merely the product of two one-particle thermal weight func- 
tions (as it would be the case for two distinguishable particles), but contains 
in addition the factor {As\A^) that accounts for the quantum effects of indis- 
tinguishability. Moreover, wf'\ai^a2) cannot be written as a product of two 
functions depending only on ai and «2- respectively, since the term ( | ) 
cannot be separated in this way. This is a result of the quantum mechanical 
principle of indistinguishability. Therefore, we say that wf'^ is entangled. 

In the case of fermions, we easily find ( A_ | A_ ) = |(1 — e^l"^^"^!^). This 
expression vanishes if cti = a2- A quantum state with two identical fermions 
in the same one-particle-state is forbidden by the Pauli exclusion principle, and 
therefore does not contribute to a thermal average. In contrast, for bosons we 
have {A+ \ A+) = |(1 -|- e'l^^^^^l ), which contains the opposite sign that 
enhances the thermal weight of the quantum state with two bosons in the 
same one-particle state. 

The thermal distribution function yields the correct partition function for 
the thermal average To show this, we calculate 

J J TT TT 2 



d^Q:i ^-|ai|2(e/3'^'^-l)^ 



TT 

]2„ i2. 



1/ 1 \' 1 

+ e- 



2 Ve^^ - ly 2 Ve2/^^ - 1 
1 'zW(/3)2 + £ZW(2/3) 

which is a correct, well-known recursion relation for the two-particle partition 
function. We indicate that the second integral is solved most easily by a change 
of variables from ai, a2 to «+ = + 0:2), = '^i^i ~ 0(2), which leads 

to a separation of the double integral. 
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2.2 Modification of the equations of motion 



The Nose-Hoover equations of motion for the coherent states parameters that 
we will investigate read 

d pi d p2 

-i-ri = — , -T.r2 = — , (5) 
dt m at m 

d 2 Pvi 

d 2 Pm 

—P2 = -mu; r2 - P2-pr- ■ 
at Q2 

The time dependence of the pseudofriction coefficients has to be determined 
in a procedure analog to the case of a single particle, i. e. we require that the 
desired distribution function 

fPiai,a2,Pr,„Pr^) oc wf{ai,a2)exp (^-/^(|^ + 
= e-"(l±e-'') 

is a stationary solution of a generalized Liouville equation on the phase space 
with elements x — {ai, Oi2,Prii,Pri2) = (n^Pi? '"2,P2,Pr?i,P??2)- The abbreviations 
U and V are defined as 

f/ = (KP + KP)(.*-i) + ;3(|- + 4) (7) 

I 12 muj , ,0 1 / \2 

V^\a,-a2\ -^{ri-r2) + ^(P^ ' P2) ■ 

As in the one-particle and regarded as classical pseudofriction 

coefficients. 

For the Liouville equation, we calculate 

^/f (ai,a2,p,„p,J = -{U + V)f^ + Ve-'' (8) 
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and 



d \ f d d d d d d 

-Pm I 

2 / 

(9) 



dx J \dri dr2 dpi dp2 dpr,^ 



— I ^v'^ 
~ Qi Q2' 

using the equations of motion (5) and imposing, as common in the present 
context [6], the constraint dpr,Jdp^^ — 0. 

After further transformations, we obtain the following equations of motion 
for the pseudofriction coefficients from a comparison of the coefficients of the 
terms Pr,i/Qi and ^,,2/^2 on both sides of the Liouville equation: 



Prn-li\-^^. ^-£P2- 



(3 \m hw mhu) + si J 

These equations, along with the set of equations (5), form a genuine quan- 
tum Nose-Hoover thermostat for two identical quantum particles. The part 
— 1 , i = 1, 2 , in the equations of motion is familiar from the dy- 
namics of a single thermostatted particle. However, in the present case of two 
particles we find additional terms that reflect the effects of Bose-attraction 
and Pauli-blocking directly in the thermostatted dynamics, see next section. 
The set of equations of motion (5), (10) conserves the quantity 

The analogous approach using a KBB-scheme starts with the set of equations 



_Pi 
dt m 
—p^ = -mwVi - gp^{Ci)Gp^{ri,pi) 

d_ _p2 

dt m 
d 



ri-'±-9r,i^i)FrAri:Pi) (12) 



'^2 = - 9ni^2)Fr2{r2,P2) 



_P2 = -TTVjJ r2 - gp^{C2)Gp^{r2,P2) , 
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with the conserved quantity 



^ l^y-^ 

The functions g^.^ , Qp^ , , 9p2 are chosen such that the canonical distribution 
of the demons can be normahzed. The functions F^^, F^j, Gpj, Gp^ are ar- 
bitrary. This scheme has the obvious advantage that positions and momenta 
are treated symmetrically, i. e. pseudofriction coefficients are present in all 
equations of motion. The time dependence of the pseudofriction coefficients 
that is obtained from the Liouville equation in the phase space with elements 
X = {n,Pi,r2,P2,6,6,Ci,C2} reads 

d_. _^(Pir dGp, , Pi - P2 1 \ . 





( ~Gpi 
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dGp, 








dpi 




( ~^P2 

\m. 
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dGp2 








dp2 



d, _ P2^ e^^-l dGp, p,-p2 1 

jIS2 — —5- I ~'^P2 ~. aZ S'-fp2 



dV (3 \m ''^ hoj dp2 mhoj + el 

dt p \ njjj ori e^+£l^ 
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dri 
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dFr, 






dr2 



The effects of Pauli-blocking and Bose-attraction are now present in the equa- 
tions of motion of all pseudofriction coefficients. 



3 Results 



3.1 Bose-attraction and Pauli-blocking 



The different signs for bosons and fermions in the equations of motion of the 
pseudofriction coefficients stem from the respective two-particle wavefunction. 

We investigate the consequences for the movements of the particles. Obviously, 
the effects of Bose-attraction and Pauli-blocking will be most pronounced at 
low temperatures, when both particles tend to occupy the one-particle ground 
state and thereby get close to one another in phase space. 

We examine the most elementary example, the scheme given by the set of 
equations of motion (5) and (10). When two fermions are close in phase space, 
V — \ai — becomes very small and the factor l/(e^ — 1) gets very large, 
causing a strong acceleration of and p^j, and thereby of pi and p2 into 
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Fig. 1. Left panel: isothermal dynamics of two bosons, right panel: two fermions. The 
initial values of {ri,pi,r2,P2} are {0, —0.1,0,0.1} in both cases, the temperature is 
T = O.lfko, and Qi = Q2 = 0.5. The integration time is 8.1 periods of the harmonic 
oscillator for the bosons, and 6.5 for the fermions. The time distance between the 
symbols is 0.013. The figure illustrates the effect of the Bose attraction and the 
Pauli-blocking in the dynamics: While the bosons stay close to each other for a 
longer time, the fermions are immediately driven away from each other due to the 
exclusion principle. 

opposite directions due to the different signs. Effectively, a close approach 
of the particles in phase space, corresponding to y — > 0, is avoided by the 
dynamics. Thus, in the case of fermions, the additional terms in the equations 
of motion (10) act like a repulsive force. 

In the case of bosons, the opposite signs in equations (5), (10) cause an accel- 
eration of the parameters in the direction of one another, favoring a "meeting" 
of the particles in phase space. So the case = is not excluded at all; on 
the contrary, it is aimed at with the maximum value of l/(e^ + !)• Fig- 1 il- 
lustrates the typical behaviour of two identical particles neighboring in phase 
space. 

In essence, the resulting effect of indistinguishability on the thermostatted 
dynamics of identical quantum particles looks like an attractive or repulsive 
interaction, although we treat a system of non-interacting particles. The in- 
teraction which is of purely statistical origin is mediated by the influence of 
the pseudofrictional forces. Nevertheless, one can hope that this statistical in- 
teraction influences the ergodicity of the system. As an extreme example, one 
could think of thermalizing only one particle and hope that the second one 
thermalizes due to the interaction, see next section. 

3. 2 Ergodicity 

The classical Nose-Hoover method applied to a single particle in a harmonic 
oscillator features ergodicity problems that are encountered in the quantum 
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case, too, since the overall characteristics of the equations of motion are sim- 
ilar in both cases. However, there is no classical counterpart of the quantum 
statistical interaction between identical particles. Therefore, the question of 
whether the more complex equations of motion presented in the precedent 
section lead to improved ergodic behavior deserves attention. 

For the different methods, we look at marginal distributions of the thermal 
weight function /J^^. As an example, we give the analytical expression of the 
marginal distribution of ri. 

1 f dpi dr2dp2 (2), 



^'^''^ ^WJ^h (15) 



1 fmuj 




22|l(e/3'i'^-l)r2 

1 -Ifisinh(^?!u;)r2 



3.2.1 Nose-Hoover and Nose-Hoover chain method 

Firstly, we present results that are obtained using the equations of motion (5), 
(10), i. e. the original Nose-Hoover technique generahzed to the case of two 
identical fermions. A detailed investigation of these equations of motion for a 
wide temperature range shows that the system is not ergodic in general. We 
find non-ergodic motion at low and intermediate temperatures T < 1.2. How- 
ever, above that value, we find ergodic motion, and the respective marginal 
distributions are well matched by the histograms obtained by time averaging, 
see Fig. 2. This fact is to be contrasted sharply to the findings in the single- 
particle case, where the simple Nose-Hoover scheme produces non-ergodic mo- 
tion even at very high temperature values [2,3]. 

The different behavior of the two-particle dynamics may be attributed to the 
statistical interaction appearing in this case. The more complicated form of 
the equations of motion improves the ergodicity of the scheme substantially 
compared to the case of a single particle. 

As a result, this technique turns out to be not recommendable in the low- 
temperature regime. Therefore, we have investigated whether the problems of 
non-ergodicity can be overcome by standard methods that are known from 
classical molecular dynamics, e. g. with a chain of thermostats [5] . 

Using one additional thermostatting pseudofriction coefficient for each p^-, 
i = 1,2, the problems of ergodicity are immediately and reliably resolved. The 
marginal distributions are exactly matched by the respective histograms, and 
the system is ergodic within the whole large temperature range that has been 
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Fig. 2. Results of time averaging with the simple Nose-Hoover scheme for two 
fermions, left panel: T = 0.1, right panel: T = 1.2. In both cases, we used 
identical initial conditions, ri(0) = — r2(0) = 0.01, pi{0) = —^2(0) = 0.01, 
Pr;i(0) = —Pr)2{^) = 0.01, Qi = Q2 = 1-0. We present density plots and histograms 
for ri, pi, and p^^; the respective results for r2, P2, and look similar. 

investigated [13]. As in the classical case, it is remarkable that a solution to the 
seemingly inaccessible problem of ergodicity is so easily obtained. Moreover, 
from the point of view of computational time, it is relatively inexpensive, since 
only two more degrees of freedom are added, which is negligible especially if 
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one deals with larger systems. 

We note that it is not sufficient to couple a second thermostating pseudofric- 
tion coefficient to only one parameter, say . In this case, the marginals of 
Prj2 are not sampled correctly at temperature values T < 0.3. Even worse, if 
only pi is coupled to a thermostat, the dynamics of the parameters r2 and p2 
is not affected at all. Hence, the energy of particle 2 is an additional conserved 
quantity, leading to strongly non-ergodic motion. This illustrates the limits of 
the improved ergodicity due to the statistical interaction. 

3.2.2 KBB method with different coupling schemes 

In the original paper [6] , Kusnezov, Bulgac, and Bauer propose the following 
choice of functions 

1111 

QrAii) = ' ^'■2(^2) = -il , QpACi) = -^Ct , 9p2{C2) = ^Cl , (16) 

as a reliable scheme that provides ergodic behavior. Using their rules of thumb 
for the adaptation of the respective coupling constants, we find that this cubic 
coupling scheme leads to ergodic motion also in the present case. However, 
considering the low-temperature results for fermions presented in Fig. 3, it can 
be inferred that the marginal distributions of the pseudofriction coefficients 
^1, and equivalently, ^2 (not shown), that are linearly^ coupled to the system 
are not sampled properly in the range T < 0.1. The histograms presented 
in Fig. 3 clearly indicate that ^1 does not depart substantially from its initial 
value during the time evolution. Nevertheless, the thermal weight function w^^^ 
is sampled correctly, since all marginal distributions coincide with the exact 
result. However, it is obvious that the cubic coupling scheme is not ergodic in 
a part of the phase space at low temperatures. This applies to bosons as well. 

We propose to resolve this problem by coupling the coefficents ^1 and ^2 in 
precisely the same manner as Ci and (2, i- e. cubically. This appears to be a 
sensible improvement since the distributions of Ci and C2 are sampled correctly 
even at very low temperatures. The resulting equations of motion yield the 
outcome shown in the right hand panel of Fig. 3. All marginal distributions 
are now sampled correctly. 

^ Recall that the derivatives of the functions given in (16) appear in the equations 
of motion. 
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Fig. 3. Marginal distributions obtained at T = 0.05 with a cubic coupling scheme 
(left panel), and with the alternate scheme described in the text (right panel). The 
initial value of on the left panel is ^i(O) = —5.1, and the "freezing" of the demon 
variable is striking. Surprisingly, all other marginals are well reproduced. The right 
panel illustrates that the problem is resolved by the modified scheme. 
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3.3 Two-particle density 



As an example for the determination of the canonical average of a two-particle 
observable, we present results of time averaging for the two-particle density, 



.(2), 



(17) 



which can also be calculated analytically. In terms of the relative and center- 
of-mass variables X- = -^{xi — X2), a;+ = -^{xi + X2), the result is given 

by 



1 muj 



1 



(. { (^^^ 2\\ 

X (l+^^^P(" ■ ^^^^ 



In order to give a useful representation of our findings, we present results of 
time averages along with the respective analytical results for fixed values of 
x+ that have been chosen arbitrarily Fig. 4 shows an excellent agreement. 





Fig. 4. Results of time averages for the bosonic (left figure) and fermionic (right 
figure) two-particle density. The solid curves correspond to the respective analytical 
results as given by equation (18). In the fermionic case, the crosses are obtained 
from the simulations that gave the right hand panels of figure 3. The upper fine 
corresponds to the value x+ = 0, the lower line to the value a;+ = 0.894. 



4 Summary and outlook 



This article presents an extension of the powerful methods of heat bath coup- 
ling in classical MD simulations to a quantum system of identical particles. 
Surprisingly, the classical ansatz turns out to be suitable for the sampling of 
quantum entangled distribution functions. The resulting two-particle dynam- 
ics contains additional terms that act like an attractive (bosons) or repulsive 
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(fermions) force mediated by the pseudofriction coefficients. Ergodicity prob- 
lems are reliably resolved by Nose-floovcr chains or a modification of the cubic 
coupling scheme, even at very low temperature values. An analytical treat- 
ment of the case of an A^-particle Fermi system is possible, and the resulting 
equations are given in Ref . [13] . 

The most promising prospect that our method offers lies in a combination with 
approximate quantum dynamics schemes. A variety of such schemes is avail- 
able, some of which are based on the time-dependent quantum variational 
principle [14] which allows to derive approximations to the time-dependent 
Schrodinger equation. How can we combine the thermostating method devel- 
oped in this work with, c. g., Fermionic Molecular Dynamics (FMD) [9] in 
order to obtain an isothermal dynamical scheme for a complex interacting 
fermion system? Does our generalization of the classical methods to quantum 
dynamics permit to make the power of approximate quantum MD schemes 
available for the calculation of quantum canonical averages without diagonal- 
ising the full many-body Hamiltonian? This question is very timely in view of 
recent experiments investigating the behavior of trapped Fermi gases. 

In FMD, the variational trial state is a Slater determinant of single-particle 
Gaussian wave packets parametrized by mean position, mean momentum, and 
complex width. These wave packets are frequently referred to as squeezed 
states and may be regarded as generalisations of coherent states. Therefore, an 
application of the thermostats developed in the present work to FMD appears 
feasible. 

Another idea of combining FMD with a thermostat is to cool the system of in- 
terest "sympathetically" , i. e. via an interaction between particles that are kept 
at a constant temperature by a quantum Nose-Hoover chain and the physical 
system under investigation. This corresponds precisely to the experimental 
technique of "sympathetic cooling" currently employed to investigate ultra- 
cold fermionic gases [15,16]. Since the thermalizing of the particles coupled 
to a Nose-Hoover chain would correspond to a thermalizing of non-interacting 
particles, such a method is only approximately correct since we need to employ 
an interaction to enable the sympathetic cooling. 

Despite these difficulties, an important asset of an isothermal MD scheme is 
that it can possibly provide temporal information, in particular time correla- 
tion functions. Although it is not clear to which degree the extended system 
methods realistically mimic the heat bath interaction, the underlying equa- 
tions of motion are physically reasonable. Therefore, in principle, this method 
is tailor-made to model the particle dynamics at constant temperature in a 
magnetic trap. 
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